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Abstract 

The benchmark active controls technology and 
wind tunnel test program at NASA Langley Research 
Center was started with the objective to investigate the 
nonlinear, unsteady aerodynamics and active flutter 
suppression of wings in transonic flow. The paper will 
present the flutter suppression control law design process, 
numerical nonlinear simulation and wind tunnel test 
results for the NACA 0012 benchmark active control 
wing model. The flutter suppression control law design 
processes using ( 1 ) classical, (2) linear quadratic Gaussian 
(LQG), and (3) minimax techniques are described. A 
unified general formulation and solution for the LQG and 
minimax approaches, based on the steady state differential 
game theory is presented. Design considerations for 
improving the control law robustness and digital 
implementation are outlined. It was shown that simple 
control laws when properly designed based on physical 
principles, can suppress flutter with limited control power 
even in the presence of transonic shocks and flow 
separation. In wind tunnel tests in air and heavy gas 
medium, the closed-loop flutter dynamic pressure was 
increased to the tunnel upper limit of 200 psf. The control 
law robustness and performance predictions were verified 
in highly nonlinear flow conditions, gain and phase 
perturbations, and spoiler deployment. A non-design 
plunge instability condition was also successfully 
suppressed. 

Introduction 

The benchmark active controls technology 
(BACT) and wind tunnel test program at NASA Langley 
Research Center was started with the objective to 
investigate the nonlinear, unsteady aerodynamics and 
active flutter suppression of wings in transonic flow. 
Under the initial wind tunnel test program, a NACA 0012 
airfoil rectangular wing, equipped with pressure 
transducers, active trailing edge control surface, and two 
spoilers were constructed for active flutter suppression 
tests. The model was mounted on a pitch and plunge 
apparatus in the NASA Transonic Dynamics Tunnel in 
order to test flutter suppression control laws and measure 
unsteady pressure distributions in nonlinear flows with 
oscillating shock and boundary layer separation. It was 
necessary to develop a flutter suppression system that 
would be stable under these flow uncertainties. 

* AIAA, Associate Fellow 


This paper describes flutter suppression control 
law design processes using (1) classical, (2) linear 
quadratic Gaussian and (3) minimax techniques. A unified 
general formulation for the linear quadratic Gaussian and 
minimax methods based on the steady state differential 
game theory is presented. Lessons learned in evaluating 
and improving the singular value based multi-input multi- 
output system robustness are described. Design 
considerations for digital implementation are outlined. A 
numerical simulation to study the effect of actuator dead- 
band, and application of the upper and lower spoiler for 
flutter suppression with the same digital control law, are 
also presented. 



Fig. 1 BACT model test set up in wind tunnel 

Wind Tunnel Model Description 

A perspective view of the BACT model test set 
up on the Pitch and Plunge Apparatus (PAPA) in the wind 
tunnel is shown in Fig. 1 . An expanded side view in Fig. 2 
shows the control surface and sensor locations. The rigid 
wing section has pitch and plunge degrees of freedom. 
The accelerometer sensors are located near the section 
leading edge ( zle ) and trailing edge (zte) at the section 
inboard. An identical pair of sensors is located at the 
section outboard as a spare. The partial span spoilers are 
located on the upper and lower surfaces, just ahead of the 
trailing edge control surface. Each of the control surfaces 
stretched over 30% of the span and 25% of the chord. The 
bending and torsion frequencies of the PAPA mounted 
NACA 0012 wing model were 3.3 Hz and 5.2 Hz 
respectively. 
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Fig. 2 NACA 0012 BACT wing on PAPA. 

Preliminary Analysis 

The preliminary analysis, control surface sizing, 
and flutter suppression control law design were based on 
the analytical state-space equations of motion of the 
BACT wing model. 14 These equations were developed 
analytically, using structural dynamic analysis and 
unsteady doublet lattice aerodynamics with rational 
polynomial approximations 5 . These linear state space 
equations consisted of 14 States (plunge, pitch, plunge 
rate, pitch rate, 3 aerodynamic states for plunge, 3 
aerodynamic states for pitch, 2 trailing edge flap actuator 
states, 2 Dryden gust states), 2 inputs (actuator command 
and gust input noise) and 7 outputs ( zte and zle 
acceleration, flap command, flap deflection, rate, 
acceleration, and gust velocity). This 14 th order state 
space equation was used for classical control law design 
and for performance simulation and verification purposes. 
For the optimal control law design purposes and for 
presentation of the design data in a concise form, the 14 th 
order state-space equations were reduced to 4 th state-space 
equations, using residualization and Schur's balanced 
reduction method 6,7 . First, it was reduced to an 8th order 
system using residualization technique, in which only the 
static part of all modes above 15 Hz were retained. The 
resulting 8th order system was then balanced and the four 
states of the system with largest balanced singular values 
were retained. A sample of the 4 th order model design data 
is presented in the Appendix. 

O pen-loop Responses 

The analytical open-loop flutter dynamic 
pressure in air was 128 pounds per square feet (psf) at a 
flutter frequency of 4.5 Hz. Fig. 3 shows the response of 
the wing trailing edge and leading edge accelerometers 
due to a 1 degree step input of the trailing edge control 
surface in air at 225 psf dynamic pressure. The primary 
plunge motion mixed with small pitch diverges rapidly. 
The unsteady lift forces oscillate about 8 lbs mean lift and 


diverges at the rate of 6 lbs/sec. The moment diverges at a 
rate of 1 lb/sec. 


Open loop TE Accl. 225 psf Open Loop LE Accl. 225 psf 



Fig. 3 Open loop transient responses in air at 225 psf. 
Frequency Responses 

The open loop frequency responses were studied 
using this 14 th order plant model, to select a possible 
candidate for feedback signal in the flutter suppression 
control law design. The Bode diagram of the trailing edge 
and leading edge accelerometers (zte) and (zle) and their 
difference (zte —zle) due to the trailing edge control 
surface excitation (8 te) in air at 225 psf at Mach 0.5, are 
shown in Fig. 4. The magnitude plots indicate 
predominant plunge response at 3.3 Hz excitation 
frequency. At 4.2 Hz excitation, the motion is a 
combination of pitch and plunge with pitch motion 
leading the plunge. The (zte —zle) represents a signal 
proportional to the pitch acceleration and can be 
integrated to provide a pitch-rate signal. Feedback of this 
signal with proper gain can provide maximum pitch 
damping at the flutter frequency. 



Fig. 4 The Bode diagrams of zte and zle and (zte — zle) 
due to 8te excitation in air at 225 psf, Mach 0.5. 
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Classical Control Law Design 

Based on this Bode plot, a classical flutter 
suppression scheme using pitch-rate proportional 
feedback from the zte and zle accelerometers was first 
devised by studying the Nyquist diagrams. The Nyquist 
diagram of the difference between trailing edge and 
leading edge accelerometers (zte - zle) due to the trailing 
edge control surface excitation (5te) in air at 200 psf, is 
shown in Fig. 5(a). The arrow indicates increasing 
frequency of excitation from 2 Hz to 6 Hz, with each * 
representing frequency increment of 1 radian/second. 
Since the open-loop plant had a pair of complex unstable 
poles, and the Nyquist contour did not encircle the -1 
point, the unit feedback closed-loop system would be 
unstable. However, if the (zte - zle) signal was integrated 
to provide a 90 degree phase lag and then used for 
feedback with sufficient gain, the Nyquist contour would 
rotate 90 degrees clockwise and then expand to encircle 
the —1 point to achieve stability. A washout filter of type 
s/(s + a) was also required, to remove any static bias that 
would otherwise be amplified by the integration. The 
series connection of integrator 1/s and washout filter was 
equivalent to a first order lag filter a/(s + a ), where s is 
the Laplace operator. 


zte, g 



(zte-zle) 



Fig. 5(a) Nyquist diagram of (zte - zle) due to Ste 
excitation in air at 200 psf, Mach 0.5 


Gain Selection 

Two types of lag filters, namely 5/(s + 5) and 
10/(s + 10) were examined. The latter was selected to 
achieve a higher phase margin at the plant input above the 
flutter frequency. Higher phase margin was desirable for 
two reasons 7 . First, the 25 Hz antialiasing filter and the 
1/200 seconds computational delay contribute about 20 
degrees of phase lag at the flutter frequency. Secondly, 
with increasing dynamic pressure, the actuators may have 
additional unknown phase lag, as the control surface 


moves against higher aerodynamic loads. The Nyquist 
diagram of the difference between trailing edge and 
leading edge accelerometers (zte - zle) with 10/(s + 10) 
lag filter and a gain KR = 500 due to the trailing edge 
control surface excitation (Ste) in air at 200 psf, Mach 0.5, 
is shown in Fig. 5(b). The unit circle is also shown. 
Because the Nyquist contour encircled the -1 point, the 
unit feedback closed loop system would be stable. As 
desired, the phase margin at the plant input above the 
flutter frequency was about 60 degrees, but the phase 
margin below the flutter frequency was only 20 degrees. 
Preliminary analysis indicated that this basic simple 
control law 1 

8 te = 500 — — (zte - zle) 
s+ 10 

can suppress the flutter instability in the dynamic pressure 
range from 0 to over 225 psf, both in air and in heavy gas 
medium. However, the closed loop transient responses 
and stability margins required substantial improvement. 


zte, g 



zle, g 



Fig. 5(b) Nyquist diagram of (zte - zle) with 10/(s+10) 
lag filter and a gain KR = 500, due to Ste excitation in air 
at 200 psf, Mach 0.5 . 

Root Locus 

Analysis of the root locus with pitch 
acceleration (zte — zle) feedback through a 10/(s + 10) lag 

filter with increasing gain KR = 0, 500 2500 is shown 

in Fig. 6(a). The stabilization was achieved by increasing 
the pitch model damping and lowering the plunge mode 
frequency. An additional feedback of the pitch rate 
proportional signal through a 5/(s + 5) lag filter with KR 
= 500 and increasing gain KP =0,500, .... 2500 was used 
to increase the damping and frequency separation further, 
as indicated by the root- locus diagram shown in Fig. 6(b). 
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This design strategy was equivalent to pitch-angle and 
pitch-rate proportional feedback that increased the pitch 
mode frequency and plunge mode damping. 


imaginary part, 
rad/sec 


imaginary part, 
rad/sec 



real part, rad/sec 


real part, rad/sec 


Fig. 6(a) Root-locus with (zte — zle) feedback through a 
10/(s + 10) lag filter, with increasing gain KR (at left). 
6(b) Root-locus with additional pitch-rate feedback 
through a 5/(s + 5) lag filter with KR = 500 and 
increasing gain KP. The arrows indicate increasing gain. 

Pitch and Pitch-Rate Feedback Control Law 

From the root-locus study, the feedback gains 
were selected as KR = 500 and KP = 1 . Thus, the second 
order state space equations of the initial pitch and pitch- 
rate feedback control law 2 is given by 


-10 0 

5 -5 

* 

X 


+ 


[10 -10 
0 0 


K =[500 1] ‘ . 


The control law inputs are zte and zle in g unit 
and the output 8te is in degrees. The high feedback gain 
was required because the maximum (zte - zle) signal was 
only of the order 0.1 g/deg, but it provided only 2% 
settling time in 1.5 seconds. Moreover, the high gain 
resulted in a severe sensitivity problem, with respect to 
uncertainty at individual sensor output, and plant 
perturbation as indicated by the corresponding singular 
value plots at 225 psf, shown in Fig. 7. Here G, K and A 
denote plant, controller and uncertainty block transfer 
function, respectively 8,9 . This figure indicates that the 
minimum singular value g(I+KG) is 0.3 at plant input and 
o(I+GK) is only 0.01 at plant output. This means that the 
closed-loop system has little robustness to multiplicative 
perturbation 8,9 at the plant output at this dynamic pressure. 

These singular value o plots can be related to 
multivariable gain and phase margins using the universal 
gain and phase margin diagram 8 shown in Fig. 8. For 
example, minimum singular value o(I+KG) of 0.3 is 


equivalent to +20 degrees phase and +3 dB gain margins 
at the plant inputs. The singular value l/o[K(I+GK) '] is 
close to 0.005 g/degrees near 2 Hz. This means that the 
plant has very little tolerance to an additive perturbation A 
to the plant. The complex determinant locus of (I+KG) 
and its distance from the origin is a measure of its 
closeness to singularity. 




_JL 


Freq, Hz 
-(G+A)-> 


Freq, Hz 
Det(l+KG) LAW2 225 psf 


.Al 



Fig. 7 Singular value plots for multivariable stability 
margin analysis due to multiplicative and additive 
perturbation A at the plant input and output, using 
classical control law 2, at 225 psf in air. 


O Minimium singular value 



Fig. 8 Universal gain and phase margin analysis diagram 

Final Pitch and Pitch-Rate Feedback Control Law 

This lack of robustness associated with this pitch 
and pitch rate feedback control law was alleviated by 
choosing a feedback of a proper linear combination of the 
two sensors with lower gains for KR, instead of using the 
difference (zte — zle). The linear combination of zte and 
zle, which is equivalent to feeding back both pitch 
acceration (zte — zle) and plunge acceleration (zte + zle) in 
the ratio 0.7(zte— zle) + 0.3(zte + zle), appeared to 
provide a superior control law. The final classical 
feedback control law, using this combination that is 
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equivalent to (zte - 0.4 zle) feedback, along with 
reduced gains of KR = 50 and KP = 1, was analyzed 
and implemented. The basic control is shown here in 
state-space form and is denoted by classical control law 3. 



-10 

5 


0 

-5 





K}=[5° 


Response and Robustness Analysis 

The closed-loop transient responses due to 1 
degree step deflection of 5te, in air at 225 psf, at Mach 
0.5, is shown in Fig. 9. The trailing edge control surface 
shows only 0.25 degrees overshoot with a maximum rate 
of 12 degrees /sec. The lift and moment forces indicate 
about 20% load alleviation compared to the open-loop 
initial transient values shown in Fig. 3. The singular- 
value plots for determining the multivariable stability 
margins 8,9 due to multiplicative perturbation at the plant 
input and output at 225 psf, in air, are shown in Fig. 10. 
Flere G, K and A denote plant, controller and uncertainty 
block transfer function, respectively. This figure indicates 
that the minimum singular value a(I+KG) is increased to 
0.8 at plant input and at plant output a(I+GK) is increased 
to 0.3 compared to those in Fig. 7. Using Fig. 8, the 
minimum singular value o(I+KG) of 0.8 is equivalent to 
+45 degrees phase and -5 db to 12 dB gain margins at the 
plant inputs. The singular value l/o[K(I+GK) '] is also 
increased from 0.005 g/degree to 0.04 g/degree near 
flutter frequency, thus increasing the plant's tolerance to 
additive plant perturbation. The complex determinant loci 
of (I+KG) ideally should be outside the unit circle to 
achieve +6dB gain margins and +60 degrees phase 
margins. The computational delay and antialiasing filters 
added 20 degrees phase lag to nearly attain these margins. 
The singular value plots indicate that the system is stable 
with adequate singular value based multivariable stability 
margins even at the high design dynamic pressure of 225 
psf. This pressure is 97 psf above the open-loop flutter 
dynamic pressure q flutter 128 psf, representing a 175% 
increase. 


Unified Optimal Design 

Flutter suppression control law design using an 
unified (1) linear quadratic Gaussian (LQG) and (2) 
Minimax method 10 - 11 is presented next. The Minimax 
approach is analogous to the time domain H-infinity 
design 12 and is based on the steady state differential game 
formulation. The unified formulation of these optimal 
design techniques provide a basic understanding of the 
relation between them. The derivation from basic 
principles using variational principles are provided. The 


solution only requires an eigen-solver. The corresponding 
Matlab script is presented in the Appendix. 



Time, sec 


Time, sec 


Fig. 9 Closed-loop responses: control surface deflection 
and rate, lift and pitching moment due to step input 8te in 
air using control law 3 at 225 psf, at Mach 0.5 (open loop 
q„ =128 psf). 

Aflutter ‘ ’ 


-(l+A)-G-> -G-(l+A)-> 



Fig. 10 Singular value plots for multivariable stability 
margin analysis due to perturbation A, using classical 
control law 3 at 225 psf in air. 

Unified Minimax Formulation 

Consider the state space Eqs.(l-3) representing 
the nth order plant, control input u(t), disturbance w(t), 
design output y £ / and sensor output y s , where all necessary 
rank, controllability and observability conditions are 
assumed to be satisfied. 

Plant state space equations 
dx(t)/dt = F x(t) + G u(t) + G w w(t) 
and x( 0) = xq (1) 

Design output 

3 ’d(t) = H d x(t) + E du u(t) (2) 

Sensor output 

y s ( t) = H s x( t) +E sw w(t) (3 ) 


5 

American Institute of Aeronautics and Astronautics 


State-Feedback Minimax Regulator Problem 

The Minimax problem is to determine the 
plant input lift) which would minimize the quadratic 
performance index J , and find the worst plant 
disturbance w(t) and initial condition x 0 , which would 
maximize J defined in Eq. (4), 

J = ( x t Q x x + x t Q xu u +u T Q i u)dt (4) 

subject to the constraint Eq.(l) with xq t xq = 1 and 
specified W defined by, 

W = \[(w T R w w)dt (5) 

Usually, the constant weighting matrices Q x , Q u are 
unity, and Q xu = [0] in a H-infinity exposition. These are 
included herein to derive a unified general time-domain 
formulation. The cross weighting matrix Q xu originates if 
one uses y d from Eq. (2) in the performance index J to 
replace x. Then, Q x = H/Q yd H d , Q xu = H d r Q yd E du , and 
Qu is replaced by [ Q u + E du r Q yd E du ]. The significance 
of the cross weighting matrix Q xu and how it can be 
selected for pole-placement of the state regulator will be 
shown later in the state-feedback regulator subsection. 

The minimax solution is given by the stationary 
condition of the augmented performance index J 

J = T.\ 0 (* 7 Qx x + 2* r e,„« + U T Q U U - y 2 w T R H ,w)dt 

+X T \“(Fx + Gu + G w w~^)dt + y 2 W (6) 

where, y is a scalar parameter. Using the calculus of 
variation with respect to x(t), u(t), w(t) and the vector 
Lagrange multiplier X(t), the conditions for dJ=0 are 
given by Eq.(l) and Eqs. (7) to (9). 


dX/dt = - Q x x - F T X - Q xu u 

(7) 

Qu u = ~ G 2 X — Q X u X 

(8) 

y2R w w = G \\[ X 

(9) 

Solving for u(t) and wft) from Eqs. 
them in Eqs. (1) and (7), the 
conditions for J are obtained as, 

(8-9) and substituting 
necessary stationary 

w =-G« _/ (G 7 ^ + QxJx) 

(10) 

w = y~ 2 R w - ] G w t X 

(ID 


* J {F-GQXO (-GQ;'G T + y->GJt- w 'GT)\Um 

Ug, +&,&’££) <f-gq?q t j t JU(Oj 

with x(0) = xq and X(oo) = 0 (12) 

State-Feedback Regulator 

Substituting Eft) = S(t)x(t), in Eqs. (10-12), leads 
to Eqs. (13-15). The general Riccati Equation (15) is then 
solved for the unknown n by n matrix S. 

u(t) = -Qu~ ] (G T S + Q xu t ) x(t) (13) 

w(t)= y~ 2 R w _1 G w t S x(t) (14) 

dS/dt + SF + F t S + Q x - (SG + QxufQiT'fSG + Q XU ) T 
+ S(T 2 G w R w _1 G w r )S = 0 (15) 

The positive definite symmetric solution for S is obtained 
from the (2n x n) eigenvectors of the n stable eigenvalues 
of the Flamiltonian matrix inside the square bracket [ ] in 
Eq.(12). For the steady state problem (i.e. dS/dt = 0 ), 
only the steady part of the Riccati Equation (15) is solved 
in order to obtain the symmetric positive-definite matrix 
S. If the eigenvectors are partitioned into two n x n 
matrices X and A which represent the stable subspace 
eigenvectors of x and X, then S = X~ ! A. The constant 
optimal feedback gains, C 0 and C w , and the closed-loop 
system matrix are given by, 

C 0 =-Qu~ 1 (G T S + Q xu T ) (16) 

C w = T 2 Rw~ 1 G w T S (17) 

dx/dt = [F + G W C W + GC a ]x. (18) 

Using Eqs. (1), (15) and (18), it can be shown 11 that 
optimal J and W defined in Eqs. (4) and (5) are given by, 

J = 0.5 Trace [S] (19) 

W = 0.5 Trace [C W T R W C W X]. (20) 

where X is the solution of the Lyapunov Equation (21) 

[F+G W C W + GC a ] X + X[F + G W C W + GC 0 ]^ + 
x(0)x(0) T = [0] (21) 

The worst x(0) that maximizes J is given by the 
eigenvector of the maximum eigenvalue of S. The 
standard linear quadratic regulator (LQR) solution is 
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obtained when y = oo, (i.e. C w = 0). As y 2 is decreased, 
the maximum singular value of [x T Q x 1/2 u t qJ' 2 ] due to 
the disturbance w(t) is reduced. The minimum value of y 2 
for which a stable solution of Eq. (15) exists provides the 
minimax state-feedback regulator that minimizes the 
maximum singular value of [x T Q x 1/2 u T Q „ 1/2 ]. 

The State-Estimator Equation 

The derivation of coupled state-estimator 
equations using linear quadratic minimax approach is still 
a subject of research. Here the equivalent state-space 
solutions 12 of the H-infinity problem are used. The state 
estimator gain B 0 D 0 is obtained by finding the symmetric 
positive definite solution for P from the state estimator 
Riccati Eq. (24) which is dual to the state regulator 
Riccati Eq. (15). 

B 0 =-(PH s t +R wv )R v ~ 1 (22) 

D 0 = (I - T 2 PSr ! , p(PS) < y 2 (23) 

dP/dt = PF t +FP + G w R w G w r - (PH s t + R wv ) R v ~ ! 
(PH s t + R wv ) t + P(T 2 Q X )P (24) 

where R v = E SW R W E SW T and must be positive definite and 
R W v = G w R w E s J. In Eq. (23) the spectrum p[PS] must 
me less than y 2 for D a to exist. The positive definite 
symmetric steady-state solution for P in Eq. (24) is 
obtained from the (2n x n) eigenvectors of the n stable 
eigenvalues of the estimator Hamiltonian matrix, 

T (F - RJG'Hf (- H\ K'H-, +y~ 2 H 1 t Q-{H 1 

l-G w KGl -R r;Xv -(F- KrK'H 2 ) 

(25) 

which is dual to the state regulator Hamiltonian matrix 
inside [ ] in Eq.(12). If the eigenvectors are partitioned 
into two (n x n) matrices X and A, then P = X / A. The 
state estimate vector z is given by the Eq. (26). 

dz/dt = F z + G w w +Gu + D„B 0 (H2z - }>2) (26) 

The complete duality relations between the state regulator 
problem and the state estimator problem are presented in 
Table 1. 

state regulator state estimator 
F F t 

G H/ 

G w H c i T 


Qx G w R w Gj 

Qu Rv 

S P 

Qxu Rwv 

(tf-t) (t-to) 

Table 1. Duality relation between linear quadratic state - 
regulator and state-estimator equations. 

The Controller Equation 

Substituting Eqs. (13), (14), (22) and (23) in Eq. 
(26), the state-estimation feedback controller equations 

dz/dt — [F + G\ X R\v + GC 0 + l x j z. — DoBoy $ (27) 

u = C 0 z . (28) 

are obtained. The standard LQG solution is obtained 
when y = oo. The minimum value of y 2 for which a stable 
solution exists for P in Eqs. (22-24), provides the minimax 
controller, that minimizes the maximum singular value of 
the matrix [y s T Q ys 1/2 u T Q u I/2 1, for the closed loop system. 

Unified Design Procedure 

In this design, the output y,j was chosen to be the 
linear combination of the trailing edge and leading edge 
accelerometer output (zte -0.4 zle), same as that used in 
the final classical design. One advantage of this choice 
was that the plant had no transmission zeros in the open 
right half plane. Usually in a frequency domain H-infinity 
design, the plant equations are augmented with weighting 
transfer functions to scale the problem into a unit sphere. 
In this time domain formulation, the weights were chosen 
as constants using the inverse of the anticipated 
j magnitude of the quantities being weighted to set their 
2 values. Using the plant Eqs. (1-3) at 225 psf in air at 
Mach 0.5 with G w = G, the initial controller was designed 
with a large value of y 2 . The unified design procedure 
block diagram for the minimax control laws is shown in 
Fig. 1 1 . The detailed design steps are described next. 

State-feedback Regulator Design 

Initially, the maximum output of the 
accelerometer sensors were of the order 0.1 g, (see figure 
3), and control surface maximum root-mean-square 
deflection was desired to be of the order 1 degrees. Thus, 
the initial values of the weighting matrices were chosen as 
follows: Q x = [ H/Q ys H s ] . Q ys = [100], and Q u = [1] . 
It was interesting to note that, instead of setting the cross 
weighting matrix Q X u = [0] as usual practice, the cross 
weighting matrix Q X u can be selected to place all state 
regulator poles beyond a certain distance a to the left of 
imaginary axis, by using 
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Qxu - - a G u (G u T G u ) 1 Qu 


(29) 

so that in Eq. (12), the eigenvalues of the diagonal matrix 
block F are off-set by 


GQu'Qxu T =~ «I. (30) 

The control-weighting matrix Q u was 

subsequently reduced to 0.01 after a few design cycles to 
improve the regulator performance. This process of 
reducing Q u is equivalent to the state estimator loop- 
transfer recovery (LTR) technique at the (nonminimal 
phase) plant output. 



Fig.ll. Unified minimax control law design and 
evaluation procedure block diagram. 

State-estimator Design 

The state estimator was designed as a dual to the 
state-regulator with R w = 1, each diagonal elements of 
R v = 0.01 , and R wv = [0], After a few design cycles the 
performance of the combined full-order controller was 
examined, and then R w was increased to 36. Since we 
also choose G w = G, this was equivalent to asymptotic 
state-regulator loop-transfer recovery (LTR) at the (non- 
minimal phase) plant input. 

4th Order Optimal Control law 

Subsequent solutions to the state-regulator and 
state-estimator were obtained with the same choice of 
weighting matrices but for decreasing value of y-, for 
which positive definite solutions for S and P could be 
obtained. Note that feasible solutions can be obtained for 
lower values of y^ up to y^> p (PS), below which the 

8 


disturbance authority exceeded the control authority. The 
4th-order optimal control law was designed with y 2 =50 to 
avoid large bandwidth controller which was typical for a 
design with no frequency dependent weighting. Fig. 12 
shows the key singular value plots for multivariable 
stability margin analysis due to multiplicative and 
additive perturbation A at the plant input and output, using 
this minimax optimal control law, denoted as control law 
4. The minimum singular value a(I+KG) is increased to 
0.9 at plant input and at plant output g(I+GK) is increased 
to 0.5 compared to those for control law 3 in Fig. 10. 


-(l+A)G-> -G(l+A)-> 



Fig. 12 Singular value plots for multivariable stability 
margin analysis due to perturbation A using nrinimax 
optimal control law 4, at 225 psf in air. 

Open Loop 225 psf Closed Loop, Law 2 225 psf 




Fig. 13 Open-loop and closed-loop response comparison 
using classical control law 2, law 3, and minimax optimal 
control law 4 due to step input 5te, at 225 psf, in air. 

Fig. 13 shows the comparison of the open-loop 
and closed-loop responses using initial control law 2, 
classical control law 3, and minimax control law 4 due to 
trailing edge control surface step input 8te in air at 225 
psf, at Mach 0.5. The transient responses indicate that the 
classical control law 3 provided better damping with 
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lower control surface activity although the minimax 
control law 4 provided better robustness properties. 
This is the traditional trade-off between performance and 
robustness. The classical control law 3 was implemented 
and tested in wind tunnel. These test results along with 
those using two optimal control laws designed by 
Waszak 13 are presented next. 

Flutter Suppression Test Results 

The performance and robustness of the final 
design was tested using the original full plant state-space 
equations and filters required for digital implementation. 
The 25 Hz antialiasing filters 157/(s+157) were added to 
the plant output. The washout filter 5s/(.?+5) and 
computational delay were added to the controller output 
equations. The 1/200 second computational delay was 
modeled by a ( 400-s)/(400+s ) filter. Before the wind- 
tunnel test entry, the digital implementation was also 
numerically simulated. The numerical simulation block 
diagram of the control system using the final classical 
control law 3 is shown in Figure 14. This nonlinear 
simulation also included the effects of a dead-band 
present in the electro-hydraulic actuator. Application of 
the upper and lower spoiler for transonic flutter 
suppression with the same digital control law was also 
investigated using this simulation. 

The active flutter suppression control-law using classical 
design was successfully tested in air and in heavy gas 
medium at transonic speeds up to Mach 0.95. The tests in 
air indicated an increase in the flutter instability boundary 
from the open-loop dynamic pressure of 158 psf (Mach 
0.38) to the tunnel limit of 200 psf. A summary of flutter 
suppression test results in heavy gas is shown in Fig. 15. 
The solid line indicates the experimental flutter boundary, 
with the transonic dip at Mach 0.8. The tests at Mach 0.8 
indicated an increase in the flutter stability boundary from 
the open-loop dynamic pressure of 142 psf to the tunnel 
upper limit of 200 psf. A non-design plunge instability 
condition was also successfully suppressed. Classical - 
control law 3 provided superior performance and waso 
demonstrated to be stable over a multiplicative gain< 
variation from 0.25 to 7, and phase variation from -90 to2 
60 degrees. A non-design plunge instability condition was® 
also successfully suppressed. Comparison of open-loopS 
and closed-loop root mean square (RMS) responses of[ii 
trailing edge accelerometer and control surfaces using the 
present classical control law 3 and two other control laws 
designed by Waszak 13 are shown in Figs. 16 and 17, 
respectively. These two control laws used upper and 
lower spoilers as control surface, for flutter suppression. 
Fig. 16 indicates that when the system is open loop stable, 
closing the loop actually reduces the response by 30%. 
Fig. 17 indicates that that the classical control law 3 
generally requires less control activity. All three control 


laws are comparable in performance, with control law 3 
exhibiting higher stability margins. 



Fig. 14 Numerical simulation block diagram of the control 
system implementation using the final classical control 
law 3 for flutter suppression. 


Closed Loop Test Points 



Fig. 15 Open loop flutter boundary and closed loop 
transonic flutter suppression tests in heavy gas. 
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0.03 

0.02 

0.01 

0.00 

Fig. 16 Open-loop and closed-loop RMS responses with 
classical control law 3 and two control laws using spoilers 
for transonic flutter suppression in heavy gas. 
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Conclusions 

Simple classical control laws, when properly 
designed based on physical principle, can successfully 
suppress transonic flutter and provide significant stability 
robustness in presence of shock and flow separation. 
Comparable robust optimal control laws can also be 
designed using a new generalized unified minimax 
formulation. Verification and improvement of the 
multivariable system stability robustness to unstructured 
perturbations at the plant, input and output were important 
steps in such a design process. Wind-tunnel tests in air 
and heavy gas indicated an increase in the transonic 
flutter dynamic pressure to the tunnel limit upper limit of 
200 psf. The control law robustness and performance 
predictions were verified in highly nonlinear flow 
conditions, gain and phase perturbations, and spoiler 
deployment. 

0.58 



M= 0.63 M= 0.71 M= 0.77 M= 0.83 
125psf 150psf 178psf 195psf 


Fig. 17 RMS responses of the control surfaces using the 
classical control law 3 and two control laws using spoilers 
for transonic flutter suppression in heavy gas. 
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Appendix 

Reduced 4 th order state space equations in air at 225 psf. 
and Matlab script for unified minimax formulation and 
solution 

F = [ -1.6073 21.0010 0. 0. 

-21.0010 -1.6073 0. 0. 

0. 0. 0.7515 25.1670 

0. 0. -25.1670 0.7515 ] 

[GG W ]= [ -3.8259 0.0597 

12.7130 0.2720 

-2.2202 -0.1107 

4.1351 -0.1745 ] 
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H d = [ -0.0517 -0.0132 -0.0668 0.0063 

-0.0542 -0.0090 -0.0482 -0.0016 

0. 0. 0. 0. 

0. 0. 0. 0. 

0. 0. 0. 0. 

10.3780 0.6369 7.9924 0.0671 

0.3897 -0.9373 -2.7625 0.9909 ] 

[E dll E dw ] = [0.0440 0.0002 

0.0421 0.0004 

1.0000 0. 

50.0000 0. 

0. 0.0968 

0.5758 -0.0004 

0.0358 -0.0003 ] 

% Matlab script for unified formulation and solution 
% xdot = F x + Gw w + G u 
% yd = Hd x + Edw w + Edu u Design output 

% ys = H x + Esw w + Esu u Sensor output 

% 7 design output [zte zle dte ddte gust lift moment] 

% — 

Qll = h'*qh*h+hd'*qhd*hd; 

Q22 = [q2]+[Edu'*qhd*Edu]; 

% pole placement using cross weight Q12 
alpha = 3.0; 

Q12=-alpha*g*((g'*g)\Q22)+[hd'*qhd*Edu]; 

WW = [Qll Q12 ; Q12' Q22]; 

% Generalized LQR Weights 
Ru = [gw*rw*gw'+nu*g*g']; 

Rv = [rv]+[Esw*rw*Esw']; 

Rwv = [gw*rw*Esw']; 

% generalized EST Weights 
VW = [Ru Rwv ; Rwv' Rv ]; 

% Generalized DESIGN with Edw=0 
[af,bf,cf,df]=lqg(f,g,h,Edu,WW,VW); 

% 

% STATE REGULATOR 

% Check if Qll is positive semi-definite and symmetric 
if any(eig(Qll) < -10*eps) I (norm(Ql l'-Qll,l)/norm 
(Q1 1,1) > eps) 

error('Qll must be symmetric positive semi-def), end 
% Check if Q22 is positive definite and symmetric 
if any(eig(Q22) <= 0) I (norm(Q22'-Q22,l)/norm(Q22,l) 
> eps) 

error('Q22 must be sym. positive definite'), end 

% 

% Construct Hamiltonian 

Hm=[(f-g*[Q22\Q 1 2']) -g* [Q22\g']+gw* [rw\gw']/mu 
-Q1 1-Q12*LQ22\Q12’] -(f-g* [Q22\Q 12’] )’] ; 

[v,d]=eig(Hm); 

% 

% Sort eigen vector of neg eigenvalues 
d = diag(d); 

[d, index] = sort(real(d)); 
if (~( (d(n)<0) & (d(n+l)>0) )) 


error('Can"t order eigenvalues'), end 
% select vectors with negative eigenvalues 
chi = v(l:n,index(l:n)); 
lambda = v((n+l):(2*n),index(l:n)); 

S = real(lambda/chi); 

% 

% Positive feedback gain ku 

ku=-Q22\(g'*S+Q12'); 

kw=[rw\gw']*S/mu; 

%~ 

% closed loop plant f = fcl 
fcl=f+g*ku+gw*kw; 

X=lyap(fcl,eye(n)) % assume xo*xo'=I 

Xg=lyap(fcl,gw*gw'*36) 

W=0.5*trace(kw'*rw*kw*Xg) 

Jo=0.5*trace(S) 

Uu=trace(ku'*ku*Xg) 

Y cov=Hd*Xg*Xg'*Hd' 

%- 

% H-inf ESTIMATOR DESIGN 
% define Hamiltonian Jam 
Ru = [gw*rw*gw'+nu*g*g']; 

J = [ (f-[Rwv/Rv]*h)' [-h'/Rv] *h+[hd'/qhd] *hd/mu 

-Ru-[Rwv/Rv]*Rwv' -(f-[Rwv/Rv]*h) ]; 
[v,d]=eig(J); 

% — 

% sort eigenvector of stable eigenvalues 
d = diag(d); [d,index] = sort(real(d)); 

% sort on real part of eigenvalues 
if (~( (d(n)<0) & (d(n+l)>0) )) 
error('Can"t order eigenvalues'), end 

% 

% select vectors with negative eigenvalues 

chi = v(l:n,index(l:n)); 

lambda = v((n+l):(2*n),index(l:n)); 

P = real(lambda/chi); 
ky = -(P*h'+Rwv)/Rv; 

% check positive definiteness of kd 

% Is this matrix (I - P*S/mu) nonsingular ? 

mumin=max(abs(eig(P*S))) 

if (mu > murnin), 

kd = inv(eye(n)-(P*S)/mu); , else 

error('spectrum(P*S) > mu , increase mu'), end 

% 

eve=eig( f+kd*ky *h ) ; 

% controller structure 

Ao=f+g*ku+gw*kw+kd*ky*h; 

evc=eig(Ao); 

% H-infinity controller 
% Kop = [Ao -kd*ky -ku zeros(nc,ns)]; 

[fc, gc ,hc, ec] = feedback(f, g, h, e, Ao, -kd*ky ,-ku, 
zeros(nc, ns)); 

% 
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